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Abstract Hawkes processes are a particularly interesting class of stochastic process that have been 
applied in diverse areas, from earthquake modelling to financial analysis. They are point processes 
whose defining characteristic is that they ‘self-excite’, meaning that each arrival increases the rate of 
future arrivals for some period of time. Hawkes processes are well established, particularly within the 
financial literature, yet many of the treatments are inaccessible to one not acquainted with the topic. 
This survey provides background, introduces the field and historical developments, and touches upon 
all major aspects of Hawkes processes. 


1 Introduction 

Events that are observed in time frequently cluster natually. An earthquake typically increases the 
geological tension in the region where it occurs, and aftershocks likely follow [1 . A fight between rival 
gangs can ignite a spate of criminal retaliations (2- Selling a significant quantity of a stock could 
precipitate a trading flurry or, on a larger scale, the collapse of a Wall Street investment bank could 
send shockwaves through the world’s financial centres [3]. 

The Hawkes process (HP) is a mathematical model for these ‘self-exciting’ processes, named after 
its creator Alan G. Hawkes [J. The HP is a counting process that models a sequence of ‘arrivals’ of 
some type over time, for example, earthquakes, gang violence, trade orders, or bank defaults. Each 
arrival excites the process in the sense that the chance of a subsequent arrival is increased for some 
time period after the initial arrival. As such, it is a non-Markovian extension of the Poisson process. 

Some datasets, such as the number of companies defaulting on loans each year 0, suggest that the 
underlying process is indeed self exciting. Furthermore, using the basic Poisson process to model say 
the arrival of trade orders of a stock is highly inappropriate, because participants in equity markets 
exhibit a herding behaviour, a standard example of economic reflexivity [6]. 

The process of generating, model fitting, and testing the goodness of fit of HPs is examined in 
this survey. As the HP literature in financial fields is particularly well developed, applications in these 
areas are considered chiefly here. 
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Fig. 1: An example point process realisation {ti,t 2 , • • •} and corresponding counting process N(t). 


2 Background 

Before discussing HPs, some key concepts must be elucidated. Firstly, we briefly give definitions for 
counting processes and point processes, thereby setting essential notation. Secondly, we discuss the 
lesser-known conditional intensity function and compensator, both core concepts for a clear under¬ 
standing of HPs. 


2.1 Counting and point processes 

We begin with the definition of a counting process. 

Definition 1 (Counting process) A counting process is a stochastic process (N(t) : t > 0) taking 
values in No that satisfies A(0) = 0, is almost surely (a.s.) finite, and is a right-continuous step 
function with increments of size +1. Further, denote by ('H(u) : m > 0) the history of the arrivals up 
to time u. (Strictly speaking %{•) is a filtration, that is, an increasing sequence of cr-algebras.) 

A counting process can be viewed as a cumulative count of the number of ‘arrivals’ into a system 
up to the current time. Another way to characterise such a process is to consider the sequence of 
random arrival times T = {Ti,T 2 ,... } at which the counting process N(-) has jumped. The process 
defined as these arrival times is called a point process, described in Definition [ 2 ] (adapted from [7]); 
see Fig. [l] for an example point process and its associated counting process. 

Definition 2 (Point process) If a sequence of random variables T = {Ti,T 2 ,...}, taking values in 
[0, 00 ), has P(0 < T\ < T 2 < • ■ •) = 1, and the number of points in a bounded region is a.s. finite, then 
T is a (simple) point process. 

The counting and point process terminology is often interchangeable. For example, if one refers to 
a Poisson process or a HP then the reader must infer from the context whether the counting process 
N(-) or the point process of times T is being discussed. 
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One way to characterise a particular point process is to specify the distribution function of the 
next arrival time conditional on the past. Given the history up until the last arrival u, H{u), define 
(as per 0) the conditional c.d.f. (and p.d.f.) of the next arrival time T k+1 as 

F*(t\H(u))= f P(T fc+1 e [s,s + ds] \H{u)) ds = [ f*(s\H(u)) ds. 

J U J U 

The joint p.d.f. for a realisation • • • , £&} is then, by the chain rule, 

k 

= p/*(tiiw(ti_i)). (i) 

i=i 

In the literature the notation rarely specifies 'H(-) explicitly, but rather a superscript asterisk is 
used (see for example 0)- We follow this convention and abbreviate F*(t\T-i{u)) and f*{t\H(u)) to 
F*(t) and /*(f), respectively. 

Remark 1 The function f*(t) can be used to classify certain classes of point processes. For example, 
if a point process has an f*{t) which is independent of R(t ) then the process is a renewal process. 


2.2 Conditional intensity functions 


Often it is difficult to work with the conditional arrival distribution f*(t). Instead, another characteri¬ 
sation of point processes is used: the conditional intensity function. Indeed if the conditional intensity 
function exists it uniquely characterises the finite-dimensional distributions of the point process (see 
Proposition 7.2.IV of 0). Originally this function was called the hazard function m and was defined 
as 


A *(i) 


/*(*) 

1 - F*(t) ' 


( 2 ) 


Although this definition is valid, we prefer an intuitive representation of the conditional intensity 
function as the expected rate of arrivals conditioned on T-l(t): 


Definition 3 (Conditional intensity function) Consider a counting process N(-) with associated 
histories %(•)■ If a (non-negative) function A *(t) exists such that 


A*(t) = lim 
hi o 


E[N(t + h)~ N(t) | n(t)\ 
h 


which only relies on information of N{-) in the past (that is, A*(t) is 'H(t)-measurable), then it is called 
the conditional intensity function of N(-). 


The terms ‘self-exciting’ and ‘self-regulating’ can be made precise by using the conditional intensity 
function. If an arrival causes the conditional intensity function to increase then the process is said to 
be self-exciting. This behaviour causes temporal clustering of T. In this setting A*(f) must be chosen to 
avoid explosion , where we use the standard definition of explosion as the event that N(t) — N(s) = oo 
for t — s < oo. See Fig. [2] for an example realisation of such a A *(t). 

Alternatively, if the conditional intensity function drops after an arrival the process is called self¬ 
regulating and the arrival times appear quite temporally regular. Such processes are not examined 
hereafter, though an illustrative example would be the arrival of speeding tickets to a driver over time 
(assuming each arrival causes a period of heightened caution when driving). 
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2.3 Compensators 

Frequently the integrated conditional intensity function is needed (for example, in parameter estima¬ 
tion and goodness of fit testing); it is defined as follows. 

Definition 4 (Compensator) For a counting process N(-) the non-decreasing function 

A(t) = I A*(s) ds 

Jo 

is called the compensator of the counting process. 

In fact, a compensator is usually defined more generally and exists even when A*(-) does not exist. 
Technically A(t) is the unique H(t) predictable function, with /1(0) = 0, and is non-decreasing, such 
that N(t) = M(t) + A(t) almost surely for t > 0 and where M(t) is an "H(t) local martingale, whose 
existence is guaranteed by the Doob-Meyer decomposition theorem. However, for HPs A*(-) always 
exists (in fact, as we shall see in Section [3j a HP is defined in terms of this function) and therefore 
Definition [4] is sufficient for our purposes. 


3 Literature review 

With essential background and core concepts outlined in Section [2j we now turn to discussing HPs, 
including their useful immigration-birth representation. We briefly touch on generalisations, before 
turning to a illustrative account of HPs for financial applications. 


3.1 The Hawkes process 

Point processes gained a significant amount of attention in the field of statistics during the 1950s and 
1960s. First, Cox [lOj introduced the notion of a doubly stochastic Poisson process (now called the Cox 
process) and Bartlett [ 1 HI 12 llT3| investigated statistical methods for point processes based on their 
power spectral densities. At IBM Research Laboratories, Lewis m formulated a point process model 
(for computer failure patterns) which was a step in the direction of the HP. The activity culminated in 
the significant monograph by Cox and Lewis m on time series analysis; modern researchers appreciate 
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this text as an important development of point process theory since it canvassed their wide range of 
applications [9} p. 16]. 

It was in this context that Hawkes [J set out to bring Bartlett’s spectral analysis approach to a 
new type of process: a self-exciting point process. The process Hawkes described was a one-dimensional 
point process (though originally specified for ieias opposed to t £ [0, oo)), and is defined as follows. 


Definition 5 (Hawkes process) Consider (N(t) : t > 0) a counting process, with associated history 
(fH(t) : t > 0), that satisfies 


P (N(t + h)~ N(t) = m | H{t)) 


\*(t)h + o(h), m = 1 
o(h ), m > 1 

1 — A *(t) h + o(h ), m = 0 


Suppose the process’ conditional intensity function is of the form 


A*(t) = A + 


f 


fi{t — u) dN(u) 


(3) 


for some A > 0 and p : (0, oo) —¥ [0, oo) which are called the background intensity and excitation function 
respectively. Assume that /r(-) ^ 0 to avoid the trivial case, that is, a homogeneous Poisson process. 
Such a process N(-) is a Hawkes process. 


Remark 2 The definition above has t as non-negative, however an alternative form of the HP is to 
consider arrivals for teM and set N(t) as the number of arrivals in (0,t]. Typically HP results hold 
for both definitions, though we will specify that this second t G R definition is to be used when it is 
required. 


(a) 


(b) 


o 


U 


t 



Fig. 3: (a) A typical Hawkes process realisation N(t), and its associated A * (t) in (b), both plotted 
against their expected values. 


Remark 3 In modern terminology, Definition [5] describes a linear HP—the nonlinear version is given 
later in Definition [6] Unless otherwise qualified, the HPs in this paper will refer to this linear form. 

A realisation of a HP is shown in Fig. [3] with the associated path of the conditional intensity 
process. Hawkes [16] soon extended this single point process into a collection of self- and mutually- 
exciting point processes, which we will turn to discussing after elaborating upon this one-dimensional 
process. 
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3.2 Hawkes conditional intensity function 


The form of the Hawkes conditional intensity function in <§ is consistent with the literature though 
it somewhat obscures the intuition behind it. Using (ti,£ 2 , • ■ •, £fc} to denote the observed sequence of 
past arrival times of the point process up to time t, the Hawkes conditional intensity is 

X*(t) = x+^2 Kt-ti)- 

ti<t 


The structure of this A*(-) is quite flexible and only requires specification of the background intensity 
A > 0 and the excitation function p(-). A common choice for the excitation function is one of expo¬ 
nential decay; Hawkes [2] originally used this form as it simplified his theoretical derivations [17] . In 
this case p(t) = ae^ 131 , which is parameterised by constants a,/3 > 0, and hence 

A*(t) = A + f ae - ^ t-s ^ dN(s) = A + ae - ^*-* -4 ^ . (4) 

-'- 00 tj<t 


The constants a and fd have the following interpretation: each arrival in the system instantaneously 
increases the arrival intensity by a, then over time this arrival’s influence decays at rate /3. 

Another frequent choice for p(-) is a power law function, giving 


A * (t) — X + 



k 

0 c+(t-s))P 


d N(a) = A +Y, 

U<t 


k 

{c+(t-ti))P 


with some positive scalars c, fc, and p. The power law form was popularised by the geological model 
called Omori’s law, used to predict the rate of aftershocks caused by an earthquake [18] • More com¬ 
putationally efficient than either of these excitation functions is a piecewise linear function as in 
IT9] . However, the remaining discussion will focus on the exponential form of the excitation function, 
sometimes referred to as the HP with exponentially decaying intensity. 

One can consider the impact of setting an initial condition A*(0) = Ao, perhaps in order to model 
a process from some time after it is started. In this scenario the conditional intensity process (using 
the exponential form of p(-)) satisfies the stochastic differential equation 

dX* (t) =/3(X — X*(t)) dt + adN(t), t> 0. 


Applying stochastic calculus yields the general solution of 

A*(£) = e _,3t (Ao — A) + A + f ae^ t-s ^ dN(s ), £>0, 

Jo 

which is a natural extension of 0 m- 


3.3 Immigration-birth representation 

Stability properties of the HP are often simpler to divine if it is viewed as a branching process. 
Imagine counting the population in a country where people arrive either via immigration or by birth. 
Say that the stream of immigrants to the country form a homogeneous Poisson process at rate A. 
Each individual then produces zero or more children independently of one another, and the arrival of 
births form an inhomogeneous Poisson process. 
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Fig. 4: Hawkes process represented as a collection of family trees (immigration-birth representation). 
Squares (■) indicate immigrants, circles (•) are offspring/descendants, and the crosses (x) denote 
the generated point process. 


An illustration of this interpretation can be seen in Fig. [4] In branching theory terminology, this 
immigration-birth representation describes a Galton-Watson process with a modified time dimension. 
Hawkes j21] used the representation to derive asymptotic characteristics of the process, such as the 
following result. 

Theorem 1 (Hawkes process asymptotic normality) If 

roo roo 

0 < n := / p(s) ds < 1 and / s/x(s) ds < oo 

Jo Jo 

then the number of HP arrivals in (0, t\ is asymptotically (t —> oo) normally distributed. More precisely, 
writing A(0,t] = JV(t) — 7V(0), 

. MoH-x/d-,) ,) 

( y«/(i-»)3 ~ ) 

where $>(■) is the c.d.f. of the standard normal distribution. 

Remark 4 More modern work uses the immigration—birth representation for applying Bayesian tech¬ 
niques; see, for example, [22] . 

For an individual who enters the system at time ti € R, the rate at which they produce offspring 
at future times t > ti is fi(t — ti). Say that the direct offspring of this individual comprise the first- 
generation, and their offspring comprise the second-generation, and so on; members of the union of all 
these generations are called the descendants of this ti arrival. 

Using the notation from 1231 Section 5.4], define Zi to be the random number of offspring in the ith 
generation (with Zq = 1). As the first-generation offspring arrived from a Poisson process Z\ ~ Poi(n) 
where the mean n is known as the branching ratio. This branching ratio (which can take values in 
(0, oo)) is defined in Theorem [l] and in the case of an exponentially decaying intensity is 

/•OO 

n= ae~^ s ds = S. (5) 

Jo P 

Knowledge of the branching ratio can inform development of simulation algorithms. For each 
immigrant i, the times of the first-generation offspring arrivals—conditioned on knowing the total 
number of them Z\ —are each i.i.d. with density p.{t — ti)/n. Section [6] explores HP simulation methods 
inspired by the immigration-birth representation in more detail. 
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The value of n also determines whether or not the HP explodes. To see this, let g[t) = E[A*(t)]. 
A renewal-type equation will be constructed for g and then its limiting value will be determined. 
Conditioning on the time of the first jump, 


ff (i)«E[A*(t)] = E 


A + f p(t — s) &N(s) = A + f g(t — s) E[dlV(s)]. 

Jo Jo 


In order to calculate this expected value, start with 


x * (s) = Um E[N(8 + h)-N(s)\H(a)] = E[dN( a )\H(s)] 
^ ' hio h ds 


and take expectations (and apply the tower property) 

E[E[dM(s)|W(s)]] _ E[dlV(s)] 


g(s) = E[A*(s)] = 


ds 


ds 


to see that 

Therefore 


E[dA r (s)] = <?(s) ds . 


g(t) = A+ / - s) g(s) ds = A + f g(t - s) /r(s) ds . 

Jo Jo 


This renewal-type equation (in convolution notation is g = A + g * fj) then has different solutions 
according to the value of n. Asmussen [23] splits the cases into: the defective case (n < 1), the proper 
case (n = 1), and the excessive case (n > 1). Asmussen’s Proposition 7.4 states that for the defective 
case 

g(t) = E[A*(t)] —> — , as t —> oo . (6) 

However in the excessive case, A*(t) — > oo exponentially quickly, and hence N(-) eventually explodes 

a.s. 

Explosion for n > 1 is supported by viewing the arrivals as a branching process. Since E [Zi] = n l 
(see Section 5.4 Lemma 2 of [23]), the expected number of descendants for one individual is 


E 




oo OO 

= £ E [Z l ]=Y j n 

i=l i= 1 



n < 1 
n > 1 


Therefore n > 1 means that one immigrant would generate infinitely many descendants on average. 

When n € (0,1) the branching ratio can be interpreted as a probability. It is the ratio of the 
number of descendants for one immigrant, to the size of their entire family (all descendants plus the 
original immigrant); that is 


E 


|2^i=i - 


1+EKrZj 1 + t^ ^ • 

Therefore, any HP arrival selected at random was generated endogenously (a child) with probability 
(w.p.) n or exogenously (an immigrant) w.p. 1 — n. Most properties of the HP rely on the process being 
stationary , which is another way to insist that n £ (0,1) (a rigorous definition is given in Section 3.41, 
so this is assumed hereinafter. 
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3.4 Covariance and power spectral densities 

HPs originated from the spectral analysis of general stationary point processes. The HP is stationary 
for finite values of t when it is defined as per Remark [2] so we will use this definition for the remainder 
of Subsection |3.4| Finding the power spectral density of the HP gives access to many techniques 
from the spectral analysis field; for example, model fitting can be achieved by using the observed 
periodogram of a realisation. The power spectral density is defined in terms of the covariance density. 
Once again the exposition is simplified by using the shorthand that 

d N(t) = lim N(t + h) - N(t ). 

^ 4-0 


Unfortunately the term ‘stationary’ has many different meanings in probability theory. In this 
context the HP is stationary when the jump process ( dN(t ) : t > 0)—which takes values in {0,1}—is 
weakly stationary. This means that E[dIV(t)] and Cov(dA r (t), dN(t+s)) do not depend on t. Stationarity 
in this sense does not imply stationarity of N(-) or stationarity of the inter-arrival times [25|. One 
consequence of stationarity is that A*(-) will have a long term mean (as given by §) 


A* 


K A* (i)] 


E[dA r (t)] 

dt 


A 

1 — n 


(7) 


The (auto)covariance density is defined, for r > 0, to be 


R{r) = Cov 


(dN(t) dN(t + r)\ 

y dt ’ dr J 


Due to the symmetry of covariance, R{—r) = R(t ), however R(-) cannot be extended to the whole of R 
because there is an atom at 0. For simple point processes E[(dA r (t)) 2 ] = E[diV(t)] (since d N(t) € {0,1}) 
therefore for r = 0 

E[(dN (t)) 2 ] = E[diV(t)] =Fd t. 

The complete covariance density (complete in that its domain is all of R) is defined as 

R (c) (r) = X*6(t) + R(t) (8) 


where S(-) is the Dirac delta function. 

Remark 5 Typically R( 0) is defined such that R^ (■) is everywhere continuous. Lewis [25 1 p. 357] 
states that strictly speaking 7U C ' (•) “does not have a ‘value’ at t = 0”. See |121I15| . and [4] for further 
details. 

The corresponding power spectral density function is then 


S(co) 


2n 



e- l ™R (c) (r) dr 


2n 


/ oo 

e~ iTul R{T) dr 

-OO 


(9) 


Up to now the discussion (excluding the final value of 0 ) has considered general stationary point 
processes. To apply the theory specifically to HPs we need the following result. 










10 


Patrick J. Laub et al. 


Theorem 2 (Hawkes process power spectral density) Consider a HP with an exponentially decaying 
intensity with a < /3. The intensity process then has covariance density, for r > 0, 

R(t) = ^ 2/3 ~ “) . e - tf - a)T ■ 

z(p — a) z 


Hence, its power spectral density is, Vw G M, 

A/3 


SW = 


27r(/3 — a) 


1 + 


a(2/3 — a) 

(/3 — a) 2 + uo 2 J 


Proof (Adapted from |2].) Consider the covariance density for r G R\ {0}: 


R(t) = E 

Firstly note that, via the tower property, 


d N(t) d N(t + t) 
dt dr 


- A* 


E 


d N(t) d N(t + r) 
d t dr 


= E 


= E 


= E 


E 

d N{t) 
d t 

d N(t) 
d t 


d N{t) d N(t + r) I , , 

~d t d? r ( * +T) 


E 


d N(t + r) 


dr 
A *(t + r) 


H(t + t) 


Hence (10) can be combined with © to see that R(r) equals 


E 


dt 


— A* 2 , 


which yields 


R(t) = + f h(t — v)R(v) dv 

J — OO 

pOO pT 

= X */i(r) + / p(r + v)R(v) dv + / fi(r — v)R(v) dv. 

J o Jo 


( 10 ) 


( 11 ) 


Refer to Appendix | A. 1 1 for details; this is a Wiener-Hopf-type integral equation. Taking the Laplace 

y \R(t) } (s) = _ Q A* (2/3 - a) _ 

1 2(/3 — a)(s + /3 — a) 

Refer to Appendix A.2 for details. Note that © and 0 supply A* = /3A/(/3 — a), which implies that 


transform of (111 gives 


% {R(t) \ (s) = __fO_ 

1 1 j 2(/3 -a) 2 (s + d-a) 


Therefore, 


R(t) = JT 


a/3A(2/3-a) \ _ a/3A(2/3 - a) -(/3-«)t 

2 (/3-a) 2 


2(/3 - a) 2 (s + ft - a) 
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S(uj) = 


2n 


X* + / e"‘™ R{t) dr 


1 

2tt 

1 

2tt 

J_ 

2tt 


The values of A* and Jzf (s) are then substituted into the definition given in & 

ir "-/ 

J — OO 

/>oo />oo 

A* + / e" iTU; J R(r)dr+ / e <TW H(r)dT 
JO Jo 

[a* +.£? {-R(r)} (iw) +«5f {-R(r)} (-zw)j 

a\*(2/3 — a) 


a\*(2/3 — a) 

A* + tvt-- 7T -r + 


2(/3 — a)(iw + /? — a) 2(/3 — a)(—iu> + /3 — a) 


A/3 


27r(/J — a) 


1 + 


a(2/3 — a) 
(/3 — a) 2 + ui 2 


Remark 6 The power spectral density appearing in Theorem [5] is a shifted scaled Cauchy p.d.f. 

Remark 7 As R(-) is a real-valued symmetric function, its Fourier transform S'(-) is also real-valued 
and symmetric, that is, 


s <“> = l 


A *+ I e~ iTU1 R(t) dr 


2tt 


\*+ COs(tui)R(t) dr 


and 

S+(w) := S(-w) + S(u) = 2S(u). 

It is common that 5+ (■) is plotted instead of S(-), as in Section 4.5 of [15]; this is equivalent to 
wrapping the negative frequencies over to the positive half-line. 


3.5 Generalisations 

The immigration-birth representation is useful both theoretically and practically. However it can only 
be used to describe linear HPs. Bremaud and Massoulie [26] generalised the HP to its nonlinear form: 

Definition 6 (Nonlinear Hawkes process) Consider a counting process with conditional intensity 
function of the form 

A *(t) — \P fj,(t — s) N(ds) 

where : R —> [0, oo), p : (0, oo) —> R. Then jV(-) is a nonlinear Hawkes process. Selecting = A + x 
reduces N(-) to the linear HP of Definition [ 5 ] 

Modern work on nonlinear HPs is much rarer than the original linear case (for simulation see 
pp. 96—116 of [7], and associated theory in G2). This is due to a combination of factors; firstly, the 
generalisation was introduced relatively recently, and secondly, the increased complexity frustrates 
even simple investigations. 

Now to return to the extension mentioned earlier, that of a collection of self- and mutually-exciting 
HPs. The processes being examined are collections of one-dimensional HPs which ‘excite’ themselves 
and each other. 
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Definition 7 (Mutually exciting Hawkes process) Consider a collection of m counting processes 
(Ml (■), ..., denoted N. Say {Tjj : i € {1,..., m},j £ N} are the random arrival times for each 

counting process (and l hJ for observed arrivals). If for each i = 1 ,m then JVj(-) has conditional 
intensity of the form 

m , t 

A*(f) = A i + ^2 Hjit-u) ANj(u) (12) 

for some Aj > 0 and m : (0, oo) —> [0, oo), then N is called a mutually exciting Hawkes process. 


When the excitation functions are set to be exponentially decaying, (12) can be written as 

rt 


K (t) — \ + ^ / 

3= 




-Pi,3 (.*-*) 


dNj(s) = \i + Yl Yl ai ’i e 

j = 1 t j, k ^ ^ 




(13) 


for non-negative constants {pt% j,P% j : i,j = 1 ,... , m}. 

Remark 8 There are models for HPs where the points themselves are multi-dimensional, for example, 
spatial HPs or temporo-spatial HPs [2j. One should not confuse mutually exciting HPs with these 
multi-dimensional HPs. 


3.6 Financial applications 

This section reviews primarily the work of A'it-Sahalia, et al. [28] and Filimonov and Sornette [6|. 
It assumes the reader is familiar with mathematical finance and the use of stochastic differential 
equations. 

3.6.1 Financial contagion 

We turn our attention to the moest recent applications of HPs. A major domain for self- and mutually- 
exciting processes is financial analysis. Frequently it is seen that large movements in a major stock 
market propagate in foreign markets as a process called financial contagion. Examples of this phe¬ 
nomenon are clearly visible in historical series of asset prices; Fig. [5] illustrates one such case. 

The ‘Hawkes diffusion model’ introduced by [28] is an attempt to extend previous models of stock 
prices to include financial contagion. Modern models for stock prices are typically built upon the model 
popularised by [29] where the log returns on the stock follow geometric Brownian motion. Whilst this 
seminal paper was lauded by the economics community, the model inadequately captured the ‘fat 
tails’ of the return distribution and so was not commonly used by traders [30j. Merton EH attempted 
to incorporate heavy tails by including a Poisson jump process to model booms and crashes in the 
stock returns; this model is often called Merton diffusion model. The Hawkes diffusion model extends 
this model by replacing the Poisson jump process with a mutually-exciting HP, so that crashes can 
self-excite and propagate in a market and between global markets. 

The basic Hawkes diffusion model describes the log returns of m assets (Xi(-),..., where 

each asset i = 1,..., m has associated expected return pi 6 R, constant volatility eq € R , and standard 
Brownian motion (W^(t) : t > 0). The Brownian motions have constant correlation coefficients 
{Pi,j '■ hJ = 1, • ■ ■ ,m}. Jumps are added by a self- and mutually-exciting HP (as per Definition [7] with 
some selection of constants a.,, and /?.,.) with stochastic jump sizes ( Zi(t ) : t > 0). The asset dynamics 
are then assumed to satisfy 


d Xi(t) = m dt + a t dW; x {t) + Zi(t) d Ni(t). 
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The general Hawkes diffusion model replaces the constant volatilities with stochastic volatilities 
{Vi(-), specified by the Heston model. Each asset i = has a: long-term mean 

volatility Qi > 0, rate of returning to this mean Kj > 0, volatility of the volatility v j > 0, and standard 
Brownian motion (W% (t) : t > 0). Correlation between the W? (•)’s is optional, yet the effect would 
be dominated by the jump component. Then the full dynamics are captured by 

d Xi(t) = iMdt+ VVJy) dW t X (t) + Zi{t) d Ni(t ), 


dVi(t) = Ki(0i - Vi(t )) d t + Viy/vJd dwf (f ). 

However the added realism of the Hawkes diffusion model comes at a high price. The constant 
volatility model requires 5m+3m 2 parameters to be fit (assuming Z^(-) is characterised by two parame¬ 
ters) and the stochastic volatility extension requires an extra 3m parameters (assuming Vi, j = 1,..., m 
that E[W4(-) V Wj(-) v ] = 0). In [28] hypothesis tests reject the Merton diffusion model in favour of the 
Hawkes diffusion model, however there are no tests for overfitting the data (for example, Akaike or 
Bayesian information criterion comparisons). Remember that John Von Neumann (reputedly) claimed 
that “with four parameters I can fit an elephant” [32] . 

For computational necessity the authors made a number of simplifying assumptions to reduce 
the number of parameters to fit (such as the background intensity of crashes is the same for all 
markets). Even so, the Hawkes diffusion model was only able to be fitted for pairs of markets (m = 2) 
instead of for the globe as a whole. Since the model was calibrated to daily returns of market indices, 
historical data was easily available (for exmaple, from Google or Yahoo! finance); care had to be 
taken to convert timezones and handle the different market opening and closing times. The parameter 
estimation method used by [28] was the generalised method of moments, however the theoretical 
moments derived satisfy long and convoluted equations. 


3.6.2 Mid-price changes and high-frequency trading 

A simpler system to model is a single stock’s price over time, though there are many different prices 
to consider. For each stock one could use: the last transaction price, the best ask price, the best bid 
price, or the mid-price (defined as the average of best ask and best bid prices). The last transaction 
price includes inherent microstructure noise (for example, the bid-ask bounce), and the best ask and 
bid prices fail to represent the actions of both buyers and sellers in the market. 

Filimonov and Sornette |B] model the mid-price changes over time as a HP. In particular they look 
at long-term trends of the (estimated) branching ratio. In this context, n represents the proportion 
of price moves that are not due to external market information but simply reactions to other market 
participants. This ratio can be seen as the quantification of the principle of economic reflexivity. The 
authors conclude that the branching ratio has increased dramatically from 30% in 1998 to 70% in 
2007. 

Later that year [33] critiqued the test procedure used in this analysis. Filimonov and Sornette 
[6] had worked with a dataset with timestamps accurate to a second, and this often led to multiple 
arrivals nominally at the same time (which is an impossible event for simple point processes). Fake 
precision was achieved by adding Unif(0,1) random fractions of seconds to all timestamps, a technique 
also used by m■ Lorenzen found that this method added an element of smoothing to the data which 
gave it a better fit to the model than the actual millisecond precision data. The randomisation also 
introduced bias to the HP parameter estimates, particularly of a and /?. Lorenzen formed a crude 
measure of high-frequency trading activity leading to an interesting correlation between this activity 
and n over the observed period. 
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Fig. 5: Example of mutual excitation in global markets. This figure plots the cascade of declines in 
international equity markets experienced between October 3, 2008 and October 10, 2008 in the US; 
Latin America (LA); UK; Developed European countries (EU); and Developed countries in the Pacific. 
Data are hourly. The first observation of each price index series is normalised to 100 and the following 
observations are normalised by the same factor. Source: MSCI MXRT international equity indices on 
Bloomberg (reproduced from [28]). 


Remark 9 Fortunately we have received comments from referees suggesting other very important works 
to consider, which we will briefly list here. The importance of Bowsher [33] is highlighted, as is the 
series by Chavez-Demoulin et al. [351136] . They point to the book by McNeil et al. m where a section 
is devoted to HP applications, and stress the relevance of the Parisian school in applying HPs to 
microstructure modelling, for example, the paper by Bacry et al. [38] . 


4 Parameter estimation 

This section investigates the problem of generating parameters estimates 6 = (A, 2, /?) given some 
finite set of arrival times t = {ti,t 2 , ■ ■ ■ , presumed to be from a HP. For brevity, the notation here 
will omit the 0 and t arguments from functions: L = L(6; t), l = 1(6; t), A *(t) = A*(i; i,0), and 
A(t) = A(t; t, 6) . The estimators are tested over simulated data, for the sake of simplicity and lack 
of relevant data. Unfortunately this method bypasses the many significant challenges raised by real 
datasets, challenges that caused [33] to state that 
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“Our overall conclusion is that calibrating the Hawkes process is akin to an excursion within a 
minefield that requires expert and careful testing before any conclusive step can be taken.” 

The method considered is maximum likelihood estimation, which begins by finding the likelihood 
function , and estimates the model parameters as the inputs which maximise this function. 


4.1 Likelihood function derivation 

Daley and Vere-Jones [9], Proposition 7.2.Ill] give the following result. 

Theorem 3 (Hawkes process likelihood) Let N(-) be a regular point process on [0,T] for some finite 
positive T, and let ti,...,t k denote a realisation of N(-) over [0,T]. Then, the likelihood L of N(-) is 
expressible in the form 

k r T 


L = 


[fiA-fo)] <*?(-/ A*(u)du) 


Proof First assume that the process is observed up to the time of the fcth arrival. The joint density 
function from 0 is 

k 

L = f(t!,t 2 ,...,t k ) = Y[f*(ti). 

i=1 

This function can be written in terms of the conditional intensity function. Rearrange 0 to find f*(t) 
in terms of A *(t) (as per [3D]): 

um _ /*(*) = j t F*{t) = dlog(l — F*(t)) 

A W 1 -F*(t) 1 -F*{t) 


dt 


Integrate both sides over the interval (t k ,t)\ 

r* 

>tk 


[ A* (u) du = log(l - F* (t)) - log(l -F*(t k )). 

Jtk 


The HP is a simple point process, meaning that multiple arrivals cannot occur at the same time. 
Hence F*(t k ) = 0 as T k+1 > t k , and so 


f A*(«) du = log(l — F* (t)). 
Jtk 


Further rearranging yields 

F* (t) = 1 — exp ^ — J A* (m) dw^, f* (t) = X* ( t ) exp ( — A* (w) du^ . 
Thus the likelihood becomes 

k k 


(14) 


(15) 


L = 


TT /*(ij) = TT A*(■ ti ) exp ( - f X* (w) du) = \ TT A*(tj) 1 exp T - [ X*(u ) du 

i= 1 i=L ^ J L l= l J V JO 


(16) 
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Now suppose that the process is observed over some time period [0,T] D [0,ffc]. The likelihood will 
then include the probability of seeing no arrivals in the time interval ( t k ,T]\ 


l = [n/*(t ; )](i-F*(T)). 


Using the formulation of F*{t) from (15 I, then 


L =*' [U ex p ( - f A*(w) du^ . 

„•_i JO 


The completes the proof. 


4.2 Simplifications for exponential decay 


With the likelihood function from (16), the log-likelihood for the interval [0, tfc] can be derived as 

k r t k k 


i = V)log(A *(ti))~ [ A*(m) du = ^ log(A* (ti)) — A(t k ). 
Jo 


Note that the integral over [0,ifc] can be broken up into the segments [0,ti], (ti,^], 
therefore 

A(t k ) = f A*(u)du = f A*(u)du+^^ f X*(u)du. 

Jo Jo Jti 


This can be simplified in the case where A*(-) decays exponentially: 

A(t k ) = f Adu-|-^^ f A + qe~ /:i( - M ~ tj ^ du 

*'° i= 1 ti tj<u 

k-1 r ti+1 * 


" | rt-i+ 1 

= A t k + a^2 E ( 

i =1 Jti 3 = 1 

k-1 i »ti. 




du 


*■__ r l i+i 

mi 


— A tk + OL 

i=l3=1 
k— 1 i 

i=i j=i 


,-P (u-tj) 


du 


e ~P(.U+i-tj) _ e -P(U-tj) 


Finally, many of the terms of this double summation cancel out leaving 


A{t k ) = \t k - - 


k-1 


(17) 

(tk-i,tk\, and 


-P(t k -u) _ „-P(u-ti) \ _ 




(18) 


Note that here the Hnal summand is unnecessary, though it is often included, see [33]. Substituting 


A*(-) and d(-) into (17) gives 
k 


i = log [ A+ a E e 

i =1 3=1 


-P(ti-tj) 


. OL 

- xtk + ft E 


-p(t k -ti) 


-1 . 


(19) 
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This direct approach is computationally infeasible as the first term’s double summation entails 
0(k 2 ) complexity. Fortunately the similar structure of the inner summations allows l to be computed 
with 0{k) complexity [4T1I921 . For i € {2,..., fc}, let A[i) = Xlj=i so that 


A(i) = e 




E e 

3 = 1 


i-U-i) + ^ e-Pfc-i-hA = e-PVi-ti- i)(i + A(i - 1)). 


j =i 


With the added base case of A(l) = 0, l can be rewritten as 

k k 

l = log(A + aA(i)) — Atfc + — 


>- P(tk-U ) 


-1 ■ 


( 20 ) 


( 21 ) 


Ozaki [8] also gives the partial derivatives and the Hessian for this log-likelihood function. Of 
particular note is that each derivative calculation can be achieved in order Oik) complexity when a 
recursive approach (similar to (20)) is taken [43] . 


Remark 10 The recursion implies that the joint process (N(t),\*(t)) is Markovian (see Remark 1.22 
of gi]). 


4.3 Discussion 


Understanding of the maximum likelihood estimation method for the FIP has changed significantly 
over time. The general form of the log-likelihood function (17) was known by Rubin 05]. It was applied 
to the HP by Ozaki [8] who derived (19) and the improved recursive form (21). Ozaki also found (as 
noted earlier) an efficient method for calculating the derivatives and the Hessian matrix. Consistency, 
asymptotic normality and efficiency of the estimator were proved by Ogata su¬ 


it is clear that the maximum likelihood estimation will usually be very effective for model fit¬ 
ting. However, [6] found that, for small samples, the estimator produces significant bias, encounters 
many local optima, and is highly sensitive to the selection of excitation function. Additionally, the 
0{k ) complexity can render the method useless when samples become large; remember that any it¬ 
erative optimisation routine would calculate the likelihood function perhaps thousands of times. The 
R ‘hawkes’ package thus implements this routine in C+-1- in an attempt to mitigate the performance 
issues. 

This ‘performance bottleneck’ is largely the cause of the latest trend of using the generalised 
method of moments to perform parameter estimation. Da Fonseca and Zaatour [20] state that the 
procedure is “instantaneous” on their test sets. The method uses sample moments and the sample 
autocorrelation function which are smoothed via a (rather arbitrary) user-selected procedure. 


5 Goodness of fit 

This section outlines approaches to determining the appropriateness of a HPs model for point data, 
which is a critical link in their application. 
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5.1 Transformation to a Poisson process 

Assessing the goodness of fit for some point data to a Hawkes model is an important practical consid¬ 
eration. In performing this assessment the point process’ compensator is essential, as is the random 
time change theorem (here adapted from 06]): 

Theorem 4 (Random time change theorem) Say {t\,t 2 , ■ ■ ■ , f&} is a realisation over time [0,T] from 
a point process with conditional intensity function A*(-). If A*(-) is positive over [0, T] and A(T) < oo a.s. 
then the transformed points {A(t\), A(t 2 ), ■ ■ ■, A(t]f)} form a Poisson process with unit rate. 

The random time change theorem is fundamental to the model fitting procedure called (point 
process) residual analysis. Original work m on residual analysis goes back to 08], 09j, and [50]. Daley 
and Vere-Jones’s Proposition 7.4.IV ]9j rewords and extends the theorem as follows. 

Theorem 5 (Residual analysis) Consider an unbounded, increasing sequence of time points {t\,t 2 , • ■ • } 
in the half-line (0,oo), and a monotonic, continuous compensator A(-) such that limt-xx, A(t) = oo a.s. 
The transformed sequence {t*,t 2 , ■ ■ ■} = (A(ti), A(t 2 ), ■ ■ ■ }, whose counting process is denoted N*(t), is 
a realisation of a unit rate Poisson process if and only if the original sequence {t\,t 2 , ■ ■ ■ } is a realisation 
from the point process defined by A(-). 

Hence, equipped with a closed form of the compensator from ( |18| ), the quality of the statistical 
inference can be ascertained using standard fitness tests for Poisson processes. Fig.[6]shows a realisation 
of a HP and the corresponding transformed process. In Fig.[6]vl(t) appears identical to N(t). They are 
actually slightly different (vl(-) is continuous) however the similarity is expected due to Doob-Meyer 
decomposition of the compensator. 


5.2 Tests for Poisson process 
5.2.1 Basic tests 


There are many procedures for testing whether a series of points form a Poisson process (see [15] for 
an extensive treatment). As a first test, one can run a hypothesis test to check ]T7 ~ Poi(t). 

If this initial test succeeds then the interarrival times, 




should be tested to ensure t/' : ^'Exp( 1). A qualitative approach is to create 
plot for Ti using the exponential distribution (see for example Fig. 7a). 
alternative is to run Kolmogorov-Smirnov (or perhaps Anderson-Darling) 


a quantile-quantile (Q-Q) 
Otherwise a quantitative 
tests. 


5.2.2 Test for independence 

The next test, after confirming there is reason to believe that the r t are exponentially distributed, 
is to check their independence. This can be done by looking for autocorrelation in the r, sequence. 
Obviously zero autocorrelation does not imply independence, but a non-zero amount would certainly 
imply a non-Poisson model. A visual examination can be conducted by plotting the points (V,;+i, Ui). 
If there are noticeable patterns then the r, are autocorrelated. Otherwise the points should look evenly 
scattered; see for example Fig. |7b| Quantitative extensions exist; for example see Section 3.3.3 of EH, 
or serial correlation tests in |3z]. 
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Fig. 6: An example of using the random time change theorem to transform a Hawkes process into a 
unit rate Poisson process, (a) A Hawkes process N(t) with (A ,a,0) = (0.5, 2, 2.1), with the associated 
(b) conditional intensity function and (c) compensator, (d) The transformed process N*(t), where 
t* = A(ti). 


5.2.3 Lewis test 

A statistical test with more power is the Lewis test as described by [53]. Firstly, it relies on the fact 
that if ■ ■ ■ ,t%f} are arrival times for a unit rate Poisson process then {t\/t* N ,t 2 /t * N ,..., t* N _ 1 /t* N } 

are distributed as the order statistics of a uniform [0,1] random sample. This observation is called 
conditional uniformity, and forms the basis for a test itself. Lewis’ test relies on applying Durbin’s 
modification (introduced in [53] with a widely applicable treatment by [55]). 


5.2.4 Brownian motion approximation test 

An approximate test for Poissonity can be constructed by using the Brownian motion approximation 
to the Poisson process. This is to say, the observed times are transformed to be (approximately) 
Brownian motion, and then known properties of Brownian motion sample paths can be used to accept 
or reject the original sample. 

The motivation for this line of enquiry comes from Algorithm 7.4.V of [9], which is described as 
an “approximate Kolmogorov-Smirnov-type test”. Unfortunately, a typographical error causes the 
algorithm (as printed) to produce incorrect answers for various significance levels. An alternative test 
based on the Brownian motion approximation is proposed here. 

Say that N(t) is a Poisson process of rate T. Define M(t) = ( N(t) — tT)/\/T for t G [0,1]. Donsker’s 
invariance principle implies that, as T —> 00 , (M(t) : t G [0,1]) converges in distribution to standard 
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Fig. 7: (a) Q-Q testing for i.i.d. Exp(l) interarrival times, (b) A qualitative autocorrelation test. The 
U k values are defined as U k = F(t k — t * k _ t ) = 1 — e _( - tfe_tfe - 1 ). 


Brownian motion ( B{t ) : t £ [0,1])- Fig. [8] shows example realisations of M(t) for various T that, at 
least qualitatively, are reasonable approximations to standard Brownian motion. 

An alternative test is to utilise the first arcsine law for Brownian motion, which states that the 
random time M* € [0,1], given by 

M* = arg max B(s ), 

»e[o,i] 

is arcsine distributed (that is, M* ~ Beta(l/2,1/2)). Therefore the test takes a sequence of arrivals 
observed over [0, T] and: 

1. transforms the arrivals to {t*/T, t^/T ,..., t k /T} which should be a Poisson process with rate T 
over [0,1], 

2. constructs the Brownian motion approximation M(t) as above, finds the maximiser M*, and 

3. accepts the ‘unit-rate Poisson process’ hypothesis if M* lies within the (a/2,1 — a/2) quantiles of 
the Beta(l/2,1/2) distribution; otherwise it is rejected. 

As a final note, many other tests can be performed based on other properties of Brownian motion. 
For example, the test could be based simply on noting that M( 1) ~ N(0,1), and thus accepts if 
M( 1) € \Z a / 2 , -Zi—a/ 2 ] and rejects otherwise. 


6 Simulation methods 

Simulation is an increasingly indispensable tool in probability modelling. Here we give details of three 
fundamental approaches to producing realisations of HPs. 


6.1 Transformation methods 


For general point processes a simulation algorithm is suggested by the converse of the random time 
change theorem (given in Section 5.11. In essence, a unit rate Poisson process (fj, t\, ■ ■ ■ } is transformed 
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(a) 




(c) 






(b) 




(d) 


eq 




Fig. 8: Realisations of Poisson process approximations to Brownian motion. Plots (a)-(c) use the 
observed windows of T =10, 100, and 10,000, respectively. Plot (d) is a direct simulation of Brownian 
motion for comparison. 


by the inverse compensator T(-) 1 into any general point process defined by that compensator. The 
method, sometimes called the inverse compensator method, iteratively solves the equations 


tt — [ A*(s) ds, t* k+1 -t* k = [ A*(s) ds 

Jo J tt. 


for {t\,t 2 , ■ • • the desired point process (see [56, and Algorithm 7.4.Ill of jjjj]). 

For HPs the algorithm was first suggested by Ozaki 0, but did not state explicitly any relation 
to time changes. It instead focused on (|14[), 



d u = log(1 - F*(t )), 


which relates the conditional c.d.f. of the next arrival to the previous history of arrivals {ti,t 2 , • ■ • ,t k } 
and the specified A*(t). This relation means the next arrival time T k+l can easily be generated by the 
inverse transform method, that is draw U ~ Unif [0,1] then t k +i is found by solving 



log (U). 


( 22 ) 


For an exponentially decaying intensity the equation becomes 


k 

log(£/) + A (t k+1 - t k ) - | 

^ i=1 


-E 




= 0 . 


Solving for t k +i can be achieved in linear time using the recursion of (20). However if a different 


excitation function is used then (22 I must be solved numerically, for example using Newton’s method 


], which entails a significant computational effort. 
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6.2 Ogata’s modified thinning algorithm 

HP generation is a similar problem to inhomogeneous Poisson process generation. The standard way to 
generate a inhomogeneous Poisson process driven by intensity function A(-) is via thinning. Formally 
the process is described by Algorithm 0EZJ. The intuition is to generate a ‘faster’ homogeneous 
Poisson process, and remove points probabilistically so that the remaining points satisfy the time- 
varying intensity A(-). The first process’ rate M cannot be less than A(-) over [0,T]. 

A similar approach can be used for the HP, called Ogata’s modified thinning algorithm |431I44| . The 
conditional intensity A*(-) does not have an a.s. asymptotic upper bound, however it is common for 
the intensity to be non-increasing in periods without any arrivals. This implies that for t G (Tj,Tj_|_i], 
A*(t) < A*(Tv“) (that is, the time just after Ti, when that arrival has been registered). So the M 
value can be updated during each simulation. Algorithm [2] describes the process and Fig. [9] shows an 
example of each thinning procedure. 


Algorithm 1 Generate an inhomogeneous Poisson process by thinning. 

1: procedure PoissonByThinning(T, A(-), M) 

2: require: A(-) < M on [0, T] 

3: P <-[],*<- 0. 

4: while t <T do 

5: E G- Exp(M). 

6: t G- t + E. 

7: U G- Unif(0,M). 

8: if t < T and U < A(f) then 

9: P«-[P, t]. 

10: end if 

11: end while 

12: return P 

13: end procedure 


Algorithm 2 Generate a Hawkes process by thinning. 

1: procedure HawkesByThinning(T, A*(-)) 

2: require: A*(-) non-increasing in periods of no arrivals. 

3: EG 10 -10 (some tiny value > 0). 

4: Pg-Q, i«-0. 

5: while t < T do 

6: Find new upper bound: 

7: M <— A*(f + e). 

8: Generate next candidate point: 

9: E G- Exp(M), t <— t + E. 

10: Keep it with some probability: 

11: U G- Unif(0,M). 

12: if t < T and U < A*(t) then 

13: P^[P,t], 

14: end if 

15: end while 

16: return P 

17: end procedure 












Hawkes Processes 


23 


(a) 
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Fig. 9: Processes generated by thinning, (a) A Poisson process with intensity A (t) = 2+sin(t), bounded 
above by M = 4. (b) A Hawkes process with (A, a,/3) = (1,1,1.1). Each (t,U) point describes a 
suggested arrival at time t whose U value is given in Algorithm [l] and Algorithm[2] Plus signs indicate 
rejected points, circles accepted, and green squares the resulting point processes. 


6.3 Superposition of Poisson processes 

The immigration-birth representation gives rise to a simple simulation procedure: generate the immi¬ 
grant arrivals, then generate the descendants for each immigrant. Algorithm[3]describes the procedure 
in full, with Fig. |10| showing an example realisation. 

Immigrants form a homogeneous Poisson process of rate A, so over an interval [0, T] the number of 
immigrants is Poi(AT) distributed. Conditional on knowing that there are k immigrants, their arrival 
times Ci, C 2 , •.., Cj, are distributed as the order statistics of i.i.d. Unif[0, T\ random variables. 

Each immigrant’s descendants form an inhomogeneous Poisson process. The ?'th immigrant’s de¬ 
scendants arrive with intensity p(t — C,) for t > Ci. Denote Di to be the number of descendants of 
immigrant i, then E[Z)j] = p,(s) ds = n, and hence D, '■AlL' Poi(n). Say that the descendants of the 

ith immigrant arrive at times (Ci + Ei,Ci + E 2 , ■. ■, Cj + Ej^.). Conditional on knowing Dj, the Ej 
are i.i.d. random variables distributed with p.d.f. fi(-)/n. For exponentially decaying intensities, this 
simplifies to Ej Exp(/3). 


6.4 Other methods 

This section’s contents are by no means a complete compilation of simulation techniques available for 
HPs. Dassios and Zhao [58] and Mpller and Rasmussen [59] give alternatives to the methods listed 
above. Also not discussed is the problem of simulating mutually-exciting HPs, however there are many 
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Algorithm 3 Generate a Hawkes process by clusters. 

1: procedure HawkesByClusters(T, A, a, /3 ) 

2: P^{}. 

3: Immigrants: 

4: k <— Poi(AT) 

5: Ci,C 2 ,...,C k Unif(0,T). 

6: Descendants: 

7: Di,D 2 ,...,D k ^ Poi(a//3). 

8 : for i <— 1 to k do 

9: if Di > 0 then 

10: Ei, E 2 ,.. ■, E d . 4 ^' Exp(/3). 

11: P t— P U { Ci + E \,..., Ci + 

12: end if 

13: end for 

14: Remove descendants outside [0, T]: 

15: P *— {Pi : Pi £ P, Pi < 7}■ 

16: Add in immigrants and sort: 

17: Sort(PU{Ci,C 2 ,...,C fe }). 

18: return P 

19: end procedure 



Fig. 10: A Hawkes Poisson process generated by clusters. Plot (a) shows the points generated by 
the immigrant-birth representation; it can be seen as a sequence of vertically stacked ‘family trees’. 
The immigrant points are plotted as squares, following circles of the same height and color are its 
offspring. The intensity function, with (A, a,/3) = (1,2,1.2), is plotted in (b). The resulting Hawkes 
process arrivals are drawn as crosses on the axis. 


free software packages that provide this functionality. Fig. [Tl] shows an example realisation generated 
using the R package ‘hawkes’ (see also Roger D. Peng’s related R package ‘ptproc’). 
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(a) 


(b) 



t 


Fig. 11: A pair of mutually exciting Hawkes processes, (a) The two counting processes N\ (t) and /^(t) 
with parameters: Ai = A 2 = 1, api = 01,2 = 02,1 = 02,2 = 2, /lip = P %,2 = @ 2,1 = P 2,2 = 8. (b) The 
processes’ realised intensities (note that A](t) = A 2 (f) so only one is plotted). 


7 Conclusion 


HPs are fundamentally fascinating models of reality. Many of the standard probability models are 
Markovian and hence disregard the history of the process. The HP is structured around the premise 
that the history matters, which partly explains why they appear in such a broad range of applications. 

If the exponentially decaying intensity can be utilised, then the joint process (A r (-), A*(-)) satisfies 
the Markov condition, and both processes exhibit amazing analytical tractability. Explosion is avoided 
by ensuring that a < p. The covariance density is a simple symmetric scaled exponential curve, and the 
power spectral density is a shifted scaled Cauchy p.d.f. The likelihood function and the compensator 
are elegant, and efficient to calculate using recursive structures. Exact simulation algorithms can 
generate this type of HP with optimal efficiency. Many aspects of the HP remain obtainable with any 
selection of excitation function; for example, the random time change theorem completely solves the 
problem of testing the goodness of a model’s fit. 

The use of HPs in finance appears itself to have been a self-exciting process. Ai't-Sahalia et al. 
[28], Filimonov and Sornette 0, and Da Fonseca and Zaatour [20] formed the primary sources for the 
financial part of Section [3j these papers are surprisingly recent (given the fact that the model was 
introduced in 1971) and are representative of a current surge in HP research. 


A Additional proof details 


In this appendix, we collect additional detail elided from the proof of Theorem [2] 
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A.l Supplementary to Theorem [2] (part one) 


R{t) = E 


= AE 


d N(t) 
d t 

dJV(i) 


U + f + \(t 

\ J — 00 


d t 

= AA* + E 


+ r — s) diV(s) 
d N(t) l rt+r 


- A* 


d t 


/ t-t -t 

fi{t 

-00 


+ r — s) dN(s) 


- X* 


dN(t) r*+r 
dt 


/ t-t-T 

-OO 


+ r — s) d 7 V(s) 


- A* 2 . 


Introduce a change of variable v = s — t and multiply by : 

d N(t) d N(t + v) 


R(t) = AA* +E 


f T rfr-v)- 

J — OO 


A * 2 = AA* + [ fi{r - u)E 
J — OO 


dN(t) dN(t + v) 


dt du 


dt du 

The expectation is (a shifted) R( c \v). Substitute that and © in: 

R(t) = AA* + f h(t — v)(R^ (v) + A* 2 ) du — A * 2 

J—00 

= AA* + f h(t — v) f A* 5 (u) + R(v) \ dv + A * 2 f h(t — v) dv — A* 

J — 00 J —00 

= AA* + A*/x(r) + f h(t — v)R(v) du + nA * 2 — A * 2 
J — 00 

= X*fi(r) + f — v)R(v) du + A* (A — (1 — n) A*). 

J — OO 


Using 0 yields 


A — (1 — n)X* = X — (1 — n) 


1 — n 


= 0 . 


R(t) = X*fi(r) + f ii(r — v)R(v) dv . 

J —OO 


A.2 Supplementary to Theorem [2] (part two) 

Split the right-hand side of the equation into three functions g\ , (72, and g% : 

r OO r>T 

R(r) = A*/i(r) -)- / p.(r + v)R(v) du + / /x(r — u)i?(t;) du . 


51O) 52(t) 

Taking the Laplace transform of each term gives 


53 ( T ) 


^{ffl(r)}(s)= /V"A*ae-^dr= -^-A* , 

/* OO /* OO 

j£? {^M} (s) = / e~ ST ae~P( r+v ^R(v) du dr 
io v/o 

rOO n oo 

= a e- pv R{v) e ~ T ( a +P) dr du 
^0 v/o 

/•oo 

- / e~P v R(v) dv 

'' Jo 


s + 0 
a 

s + 0 


&{R} (£), 


du — A* . 


( 23 ) 
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& { 23 (t)} (s) = -Sf {m(t)} (s)-Sf {R(t)} (s) = {^( T )} ( s ) • 


Therefore the Laplace transform of ( |23| ) 


-s? {* M} (-) = (a* + if {fl(r)} (/3) + Se{Bir)} (sj) . 


(24) 


Substituting s = /3 and rearranging gives that 



(25) 


So substituting the value of jSf {i?(r)} (/3) into ( |24[ ) means 
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